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This paper presents a discrete adjoint method for a broad class of time-dependent op- 
timization problems. The time-dependent adjoint equations are derived in terms of the 
discrete residual of an arbitrary finite volume scheme which approximates unsteady con- 
servation law equations. Although only the 2-D unsteady Euler equations are considered 
in the present analysis, this time-dependent adjoint method is applicable to the 3-D un- 
steady Reynolds-averaged Navier-Stokes equations with minor modifications. The discrete 
adjoint operators involving the derivatives of the discrete residual and the cost functional 
with respect to the flow variables are computed using a complex-variable approach, which 
provides discrete consistency and drastically reduces the implementation and debugging 
cycle. The implementation of the time-dependent adjoint method is validated by com- 
paring the sensitivity derivative with that obtained by forward mode differentiation. Our 
numerical results show that 0(10) optimization iterations of the steepest descent method 
are needed to reduce the objective functional by 3-6 orders of magnitude for test problems 
considered. 


I. Introduction 

Time-dependent optimization problems arise in many areas in science and engineering including various 
flow control applications such as controlling flow separation, airframe vibration, noise level, transition to 
turbulence, etc., as well as design optimization problems for essentially unsteady flows, including design and 
shape optimization of helicopter rotors, turbomachinery blades, aircraft wings, and other configurations. 
The overall complexity of these time-dependent optimization problems is much higher than that of steady- 
state aerodynamic optimization problems, which is one of the main reasons why this class of problems 
has not been used in real-life applications. Continuously expanding computer capabilities now allow more 
attention to be devoted to numerical solution of these class of time-dependent optimal control and design 
optimization problems. These problems can be considered as minimization of an appropriate cost functional 
(e.g., lift, drag, etc.). The resulting control laws or design variables are obtained by solving the corresponding 
time-dependent optimal control or design problems by using appropriate optimization algorithms. 

Among various optimization techniques available in the literature (see, e.g., f 1 * ]), the adjoint method 
has recently grown in popularity, rapidly becoming one of the most widely used techniques for solving a 
variety of steady and unsteady optimization problems. The adjoint methodology is particularly attractive 
for optimal control/design problems, which include a large number of control variables, yet relatively few 
constraints. In contrast to a classical forward mode differentiation approach which requires solutions as 
twice as many flow problems as the number of the corresponding control variables, the adjoint methodology 
has the advantage of computing the cost functional gradients at a fixed expense independent of the number 
of control/design variables. This property of the gradient methods based on the adjoint formulation make 
them very well suited for steady aerodynamic design optimization problems. 2-5 Although the adjoint-based 
methods have been successfully used for problems of optimal design within the steady-state aerodynamics, 
applications of the adjoint formulation to essentially time-dependent optimal control/design problems are 
still lacking. In, 6 the 2-D continuous time-dependent adjoint incompressible Navier-Stokes equations and 
optimality conditions have been derived. This continuous adjoint-based method has been successfully used 
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for solving the problem of boundary-layer instability suppression through wave cancellation. Nadarajah and 
Jameson 7 derived and applied the time accurate continuous and discrete adjoint equations to the shape 
optimization of an oscillating airfoil in an 2-D inviscid transonic flow. In, 8 a gradient method based on the 
discrete adjoint equations and the corresponding boundary conditions in the frequency domain has been 
developed. This approach significantly reduces the computational cost for shape optimization of a 3-D wing 
oscillating at a constant frequency. Note, however, that this technique is applicable only for periodic problems 
and its efficiency strongly depends on the number of harmonics in the time-dependent solution. Discrete 
adjoint-based methods operating directly in the time domain have been developed for the 2-D compressible 
Euler and Navier-Stokes equations in 9 and, 10,11 respectively. 

The time-dependent adjoint-based methods mentioned above can be divided into two groups. The first 
one 7-9,11 is developed for optimization problems with the design/control variables that are independent of 
time, while the second one 6,10 involves the design/control variables that depend on time. In this paper, 
we develop a general discrete adjoint-based optimization methodology which is directly applicable to both 
classes of problems. This time-dependent optimization methodology can be directly applied to solving a 
very broad spectrum of time-dependent optimal control problems, where the control variables are in general 
time-dependent (e.g., the displacement of an actuator diaphragm or the velocity distribution at the actuator 
orifice, etc.) and design optimization problems where the design variable are in general do not depend on 
time (e.g., shape of a helicopter rotor or an aircraft wing, etc.). 

The paper is organized as follows. In Section II, we present the continuous and discrete state equations. 
In Section III, the discrete time-dependent optimization problem is described. In Section IV, the general 
discrete time-dependent adjoint equations are derived. Section V discusses a technique for forming discrete 
adjoint operators by using complex variables. In Section VI, we present two test problems used for validating 
the developed time-dependent optimization methodology. We draw conclusions and present our plans for 
the future in Section VII. 


II. Governing Equations 


We consider the time-dependent, two-dimensional Euler equations describing the unsteady, inviscid com- 
pressible flow. The Euler equations written in the integral conservation law form are given by: 


dV\J 

dt 



ncff = 0 


(1) 


where n is the outward unit normal vector of the control volume with boundary T, V is the control volume, 
U is the vector of conserved variables averaged over the control volume, and F is the Cartesian inviscid flux 
vector. 

The governing equations (1) are discretized by using a node-centered finite- volume scheme, where solution 
values are stored at the mesh nodes. The control volume around each grid node is constructed by connecting 
the centroids of the primal- mesh cells with midpoints of the surrounding edges. The discretized Euler 
equations including the boundary conditions can be written as follows: 


Q n - Q"- 1 
At 


+ R(Q”) = 0, 


(2) 


where Q = VU, and R is the spatial undivided residual of the discretization which approximates the contour 
integral in Eq. (1). It should be noted that the above discrete formulation (2) is very general and can be 
applied to a broad class of time-dependent PDEs discretized using not only finite volume, but also finite 
difference and finite element schemes. The flux F in the discretized integral is approximated using Roe’s 
approximate Riemann solver 


F 


1 - 
2 . 


Fl + F r 


|A| (Ul — Uk) , 


(3) 


where F^ and F^ are the “left” and “right” normal fluxes at the edge midpoint, U l and are the 
“left” and “right” reconstructed values of the solution vector at the edge midpoint, obtained from some 
polynomial approximation defined on each control volume, | A | is the Roe averaged matrix. In Eq. (2), the 
time derivative has been approximated using the implicit first-order backward-difference (BDF-1) formula. 
Note that second-order BDF formula as well as higher order implicit Runge-Kutta methods can also be used 
in the present formulation with minor modifications. 
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In the present paper, we consider only inviscid flow problems because our primary objective is to develop 
a time-dependent adjoint-based optimization methodology that is applicable to a broad spectrum of nonlin- 
ear state equations. Generalization of this methodology to the unsteady Reynolds-averaged Navier-Stokes 
(RANS) equations coupled with either one- or two-equation turbulence model is quite straightforward. In 
this case, only the flux residual R should be changed, while the adjoint equations, which will be presented in 
Section IV, remain unchanged. Note, however, that for the RANS equations, questions related to robustness 
of the present adjoint-based methodology require special investigation, which are outside the scope of the 
present paper. 


III. Discrete Time-Dependent Optimization Problem 

We consider the following discrete time-dependent optimization problem: 

N 

min /(D), /(D) = £ /"(D)Ai, /"(D) = i£ bj (CT(D)) + F£ eg (U, D) (4) 

UGT'a n— 1 

where D is a vector of the control or design variables which in general depends on time, N is the total 
number of time steps over which the control D is active, Q is the solution of the unsteady, compressible 
Euler equations, F™ h ■ is a part of the cost functional that represents the flow control objective, and F™ g 
is a regularization term, typically some weighted norm of the control variable. Note that this setting of 
the problem remains valid if D does not depend on time. The above formulation (4) is very general and 
directly applicable for both time-dependent optimal control problems (e.g., active flow control via synthetic 
jet actuation) and aerodynamic design optimization for unsteady flows (e.g., design of a turbomachinery 
blade), and others. The set of admissible controls, V a , depends on specifics of the target physical system 
(e.g., how much suction and blowing can the actuators provide, or admissible length and thickness of a blade, 
etc.), but it should also ensure the existence of a solution of the optimization problem (4). 

To reduce the complexity of the optimization problem, without loss of generality, we assume that the 
objective functional / is a scalar quantity. In the present analysis, we consider the following discrete convex 
functional: 

Fo-b, = A (C£ - C£.„,.,) 2 + A (cs - CB,„J 2 (5) 

where C£ tar et and CS target are given time-dependent target lift and drag coefficients, respectively, which are 
integrals of the normal and tangential components of the stress tensor over the controlled boundary surface. 

The control variables D have a precise physical meaning (e.g., the the Mach number or angle of attack as a 
function of time, etc.) and should remain bounded and be continuous in time. These physical constraints are 
incorporated into the optimization problem through the regularization/penalty term in the cost functional 
Eq. (4), which limits the size of the control. The regularization/penalty term F r £ g is chosen as follows 

K g = H [D"] T D" + (D ra - D"” 1 ) T (D" - D”- 1 ) (6) 

where oq and «2 are nonnegative parameters that can be used to adjust the relative weights of the reg- 
ularization terms appearing in the functional (6). The particular form of the penalty term (6) limits not 
only the magnitude of the control, but also the rate, at which the control changes, to provide the necessary 
smoothness of the control. The presence of the second term in the cost functional can also be interpreted 
as a constraint on the maximum kinetic energy generated by the control system, which is directly related to 
the energy consumption required for its operation. 

It should be noted that the same penalty technique outlined above can be used to impose a more general 
nonlinear side constraints involving the state variable U. If the optimization problem (4) is subject to the 
side constraint <I>(U, D) < 0, then the following penalty term can be added to the objective functional to 
enforce this constraint: 

Er eg = a (max [0, < f>(U, D)]) 2 , (7) 

where cc is a positive user-defined parameter, and $ is a continuously differentiable function of its arguments. 
Note that the above penalty term is continuously differentiable and active only when the constraint is 
violated, i.e., when < f>(U, D) > 0. The above penalization guarantees that the constraint 4>(U, D) < 0 is met 
if (3 — > +oo. In practice, the parameter a can be increased during the iterative process to make sure that 
the side constraint is satisfied. 
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Currently, only matching objective functionals with the global extremum have been considered. Because, 
there are no spurious extrema in all test problems presented herein, the regularization term is set equal to 
zero. Though this simplified formulation works well for the Euler equations considered in this paper, for 
problems involving essentially nonlinear one- or two-equation turbulence models, the regularization term 
may play an important role and should be included into the optimization procedure. 


IV. Time-dependent Adjoint Formulation 


The discrete time-dependent optimization problem (4) is solved by using the method of Lagrange mul- 
tipliers which is used to enforce the governing equations and the corresponding boundary conditions (2) as 
constraints. The discrete Lagrangian functional is defined as follows: 

L(D,Q,A)= £ /"At + E [A"f ( Qn -2r + R ") At + [ Al ] T (Q 1 - Q in ) , (8) 

71=1 71=2 V 7 


where A is a vector of Lagrange multipliers or costate variables, n = 1 corresponds to the initial moment 
of time, Q m is the initial condition for the Euler equations, /" is the objective functional given by Eq. (5), 
and R n = R(Q", D) is the spatial undivided residual. Note that the first two terms in the Lagrangian 
are scaled by At, so that they approximate the corresponding time integrals in the continuous Lagrangian. 
Therefore, the discrete Lagrangian approaches its continuous counterpart as the number of time steps N 
increases. Furthermore, the scalar product of the costate vector and the vector of the governing equations 
in Eq. (8) can be interpreted as the integral over the computational domain, which again approximates the 
continuous Lagrangian. 

The sensitivity derivative is obtained by differentiating the Lagrangian with respect to D, which yields 
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Regrouping the terms, Eq. (9) can be recast as follows: 
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(10) 


For problems with a large number of control/design variables, it is desirable to avoid the calculation of 
<9Q/<9D in the optimization procedure, because it requires multiple solves of the primary problem. Taking 
into account that so far, no constraints have been imposed on the Lagrange multipliers, the 9Q/9D term 
can be eliminated from Eq. (10) by setting the second, third, and forth terms on the right hand side equal 
to zero, which results in the following adjoint equations for determining the Lagrange multipliers: 



<9R 
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-,T 


dQ 
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L N 


df N 

dQ N 


( 11 ) 
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At 
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<9R n 


0Q r 
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At 


(A 1 


A 2 ) 


df n 
dQ n ’ 
Of 1 
9Q 1 ’ 


2 < n < N - 1 


(12) 

(13) 


Equations (11) and (13) are initial and terminal conditions for the costate variables. Equations (12) represent 
a linear system of equations for the costate variables, which are solved backward in time. Once Eqs. (11-13) 
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have been solved, the vector of Lagrange multipliers A" can be used to evaluate the last two terms in Eq. 
(10). As a result, the sensitivity derivative can be calculated as follows: 


dL 

dD 


N ■ 1 rn X 


3R" 


3D 


A” At 



(14) 


where and are calculated by using Q" stored during the forward sweep in time. 

The minimum of the functional is found by using the steepest descent method in which each step of the 
optimization cycle is taken in the negative gradient direction 

D ™ +1) = D ™- t ™^’ ( 15 ) 

where r m is the step size for the m-th component of the vector D, k is the number of optimization cycle, D m 
is the ?n-th component of the vector D. The sensitivity derivative dL/dD in Eq. (15) is determined using 
Eq. (14) which requires the solution of the adjoint equations (11-13). During the solution of the adjoint 
equations that are integrated backward in time, the sensitivity derivative at each time step is computed and 
added to its value at the previous time step. At n = 1, the complete sensitivity derivative is available and 
used in Eq. (15) for calculating a new value of the control variable D (fc+1 / Then the entire optimization 
cycle is repeated until \L^ k+1 ^ — L^\ < e, where e is a given tolerance. This optimization algorithm is quite 
sensitive to the step size r and has been selected because of its simplicity. Other more efficient and robust 
gradient-based methods, such as conjugate gradient or quasi-Newton methods, can also be easily coupled 
with the time-adjoint formulation used in the present analysis. 

As has been mentioned above, the selection of the optimization step size r in Eq. (15) has a very strong 
effect on the number of optimization iterations required to reach an optimum solution. Large step sizes may 
result in instabilities in optimization iterations, whereas small step sizes provide stability, but drastically slow 
down the convergence. Therefore, the optimization step size in Eq. (15) is selected adaptively to maximize 
the convergence rate. In the present analysis, the following algorithm for choosing r is used: 

1) Each optimization cycle continues until a solution with a smaller integrated cost functional, 

N 

f n At, has been found. 

n — 1 

2) For a trial vector of r, the new vector D and the corresponding U are computed on all time levels. 

3) The n-th component of the vector r, r„, is changed if one of the following two events occurs: 

4) if the local cost functional increases, i.e., / t " ial > /", then r„ is decreased, r„ = 0.5r„. 

5) or if the local cost functional decreases slowly, /" > / t ” ial > 0.9/", 

then T n is changed depending on signs of j at the current and the previous 
optimization cycles. 

( \ current / \ previous 

) ( JB ~ ) < 0i then r„ is decreased, r„ = 0.5r„, 

7) otherwise, the signs are the same and r n is increased, r n = 1.5r n . 

Here, / t " ial and /" are the trial and current values of the objective functional at the time level n, and dL/dD n 
is the sensitivity derivative with respect to the control variable D n . 


V. Forming Discrete Adjoint Operators by Using Complex Variables 


As follows from Eqs. (11-13), the derivatives of R and / with respect to Q and D are required to 
form the adjoint equations and the sensitivity derivative. It is very difficult and time-consuming to obtain 
these derivatives by manually differentiating a CFD solver, especially if complicated turbulence and physical 
models are involved. Furthermore, any changes in the discretization of the governing equations, boundary 
conditions, objective functional, or physical models require additional coding and debugging, thus making 
the software development cycle extremely lengthy. To overcome these difficulties, we use an approach based 
on complex variables, which has successfully been applied to solving design optimization problems in [ 12, 13 ]. 
The key idea of this technique is to approximate the required real-valued derivatives by using the following 
formula proposed by Lyness 14 : 


df Im [f(x + ih)] 2 

dx ~ h 


(16) 
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where f(x) is a complex-valued function. In contrast to the finite-difference method, the above complex 
variable formula is robust for small h, while providing true second order accuracy. Another advantage of this 
approach is that no additional flow solves are required to evaluate this derivatives, because only the solution 
at the current time level Q n is needed to compute R and / and their perturbed values. The complex variable 
approach drastically reduces the implementation cycle and provides adjoint-based optimization capabilities 
for realistic physical and turbulence models. Note, however, that this approach is not without penalties in 
the CPU time and memory as compared with the handcoded Jacobians implementation because complex 
arithmetic is used. 



Figure 1. Comparison of the sensitivity derivatives computed using the finite difference and adjoint methods for the 
first test problem. 



Figure 2. Convergence history of the objective functional ^obj • 


VI. Numerical Results 

In this section, we present computational results demonstrating how the adjoint-based method performs 
for two time-dependent optimization problems, involving flow matching functionals. The key distinction 
between these two problems is that for the first one, a design/control variable is independent of time, while 
for the second problem, control variables depend on time. 
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Time 

Figure 3. Comparison of the optimal lift coefficient computed using the adjoint-based method with its initial and target 
values for the first test problem. 

The first problem is a minimization of a matching functional given by Eq. (5) with /?2 = 0 for the 
unsteady flow around a bump. The unsteadiness is introduced into the flow through the freestream Mach 
number which oscillates in time 

M(t) = M 0 + AMcos(wf), (17) 

where Mq is a mean value of the freestream Mach number, AM is a Mach number amplitude, and to is a 
frequency of Mach number oscillations. The thickness of the bump is set to be 10% of its chord length. 
The flow conditions used in this test problem are: M 0 = 2.0, to = 177t/ 9, and Tg na i = 1. For this test 
problem, the Mach number amplitude AM is used as a control variable. Note that this control variable 
is independent of time. The time-dependent target lift coefficient Cx t t (t) in the objective functional is 
calculated numerically by solving the unsteady Euler equations (2) with AM = 0.5. A solution obtained at 
AM = 0.1 is used as a starting point for the optimization procedure. The optimization is stopped when the 
absolute value of the difference between the current value of the Lagrangian and its value at the previous 
optimization iteration is less than 10 _8 . This test problem is solved on a 61 x 21 structured grid using a 
node-centered finite volume code that is first-order accurate both in time and space. At each time step, 
the nonlinear discrete flow equations are solved by using the Newton’s method. The adjoint equations are 
integrated backward in time and require the solution of the Euler equations to be known at all time steps, 
over which the optimization problem is solved. In the present implementation, the entire unsteady solution 
set is held in operating memory. For mid-size 2D problems integrated over O(10 2 ) time steps, the required 
solution history, can be stored on a hard drive, as has been reported in [ 9 ], In this case, the speed of I/O 
operation itself does not have a significant effect on the overall CPU time. Note, however, that for realistic 
3D nonperiodic problems, this approach can quickly become prohibitive in terms of the disk memory and 
the CPU time, so more efficient approaches are needed to solve this class of unsteady optimization problems. 

To evaluate the accuracy of calculation of dL/dD and to check the implementation of the adjoint solver, 
two different methods for computing the sensitivity derivative are used and compared with each other. The 
first method is based on a forward mode differentiation of the Lagrangian with respect to the control variable, 
which is implemented by using the complex variable approach (16). The second method uses the discrete 
adjoint formulation described in Section III. Note that at each optimization cycle, for the first method, the 
flow problem should be solved twice for each control variable. In contrast to the first method, the second 
approach requires a single solution of the Euler and corresponding adjoint equations per optimization cycle, 
regardless of the number of the control/design variables. For the forward mode differentiation, the complex 
step size has been chosen to be 1CU'. Figure 1 shows the difference between the sensitivity derivatives 
obtained using the finite difference and adjoint methods. As seen in the figure, the discrepancy is of the 
order of the round-off error, thus corroborating the validity of the time-dependent adjoint formulation. The 
history of convergence of the objective functional is presented in Fig. 2. The value of the objective functional 
drops by an order of magnitude every 2 optimization cycles. Note that this convergence behavior remains 
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Figure 4. Comparison of the target Mach number and the optimal Mach number computed using the adjoint-based 
method for the second test problem. 





10 20 30 40 50 60 

Number of Iterations 


Figure 5. 


A history of convergence of the adjoint-based optimization method for the second test problem. 


practically unchanged until the functional becomes smaller than the specified tolerance when the optimization 
was stopped. Only 10 optimization iterations were needed to reduce the objective functional by six orders of 
magnitude. To illustrate that the lift coefficient converges to its target value, time histories of the optimal, 
target, and initial lift coefficients are depicted in Figure 3. At the first optimization iteration, the maximum 
value of the time-dependent lift coefficient is about two times less that that of the target Cz,(t). After 10 
optimization iterations, the time history of the computed lift coefficient is practically indistinguishable from 
the target solution at all time steps. 

The second test problem is similar to the first one, but now, values of the freestream Mach number at 
each time step D = (Mi, . . . , Mn) t are used as control variables. The optimization procedure starts at 
M = 2.1 which is used as an initial guess. The objective functional for this test problem is given by 

Fobj = E ( P 7 - target ) 2 , 
r c 

where PJ 1 and P t ™ rget are computed and target time-dependent pressure profiles at the lower boundary of 
the computational domain. The target pressure distribution is calculated numerically by solving the same 
unsteady problem with M = 2 + 0.5 cos(177rf/9). The optimization is stopped when either the relative 
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Figure 6. Comparison of the optimal lift coefficient computed using the adjoint-based method with its initial and target 
values for the second test problem. 


change in the value of each control variable becomes smaller than 1CU 4 or the absolute value of the objective 
functional becomes smaller than 1CU 7 . As in the previous case, the adjoint equations are integrated backward 
in time and require the solution of the Euler equations to be known at all time steps, which is held in the 
operating memory. Figure 4 shows the optimal and target Mach number distributions in time. As seen 
in the figure, the time-dependent optimization method converges to the target solution on the entire time 
interval considered, thus validating the unsteady adjoint formulation. 

A history of convergence of the objective functional obtained with adjoint-based optimization technique 
is presented in Fig. 5. The total number of optimization cycles required for the adjoint-based optimization 
method to converge is an order of magnitude larger than that obtained for the first test problem. This is 
not surprising, because the dimensionality of the design space has also increased by an order of magnitude. 

To illustrate that the lift coefficient converges to its target value, the optimal Cl obtained with time- 
dependent optimization method as well as the initial and target lift coefficients are presented in Figure 6. 
The relative difference between the initial lift coefficient and its target value is of the order of 0(1), while 
the solution obtained with the adjoint-based method is almost indistinguishable from the target C l over the 
entire time interval considered. 


VII. Conclusions 

We have developed the general adjoint-based methodology for solving a broad spectrum of optimal flow 
control and design optimization problems. The methodology is directly applicable to both time-dependent 
optimization problems with control/design variables that are time-dependent and design optimization prob- 
lems with the control variables that do not depend on time. Nonlinear constraints on the control/design and 
state variables can be incorporated into the present formulation by introducing the penalty/regularization 
term in the cost functional. The discrete adjoint operators required for this formulation are computed by 
using the complex variable approach which is robust for very small step sizes, thus providing adjoint-based 
optimization capabilities for realistic physical models. The present adjoint-based methodology has been 
validated using two test problems involving flow matching functionals. Applications of this adjoint-based 
methodology to more realistic time-dependent design optimization problems involving moving and deforming 
grids is currently under investigation. Our future research will also focus on developing optimization and 
computational techniques for reducing the CPU and memory cost of the present time-dependent adjoint- 
based methodology. 
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